Update to http://gis.stackexchange.com/questions/90635/what-programs-would-allow-for-the-mapping-of-a-geoid-in-3d

Download the files with R:

    baseurl <- "http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/GIS/world_geoid"

    fs <- file.path(baseurl, c("n45w180.zip", "n45w135.zip", "n45w90.zip", "n45w45.zip",    "n45e00.zip", "n45e45.zip", 
    "n45e90.zip", "n45e135.zip", "n00w180.zip", "n00w135.zip", "n00w90.zip", 
    "n00w45.zip", "n00e00.zip", "n00e45.zip", "n00e90.zip", "n00e135.zip", 
    "s45w180.zip", "s45w135.zip", "s45w90.zip", "s45w45.zip", "s45e00.zip", 
    "s45e45.zip", "s45e90.zip", "s45e135.zip", "s90w180.zip", "s90w135.zip", 
    "s90w90.zip", "s90w45.zip", "s90e00.zip", "s90e45.zip", "s90e90.zip", 
     "s90e135.zip"))

    for (i in seq_along(fs)) download.file(fs[i], basename(fs)[i], mode = "wb")

     ## unzip
    for (i in seq_along(fs)) utils::unzip(basename(fs[i]))

    ## build the next string with R too
    txt <- gsub(".zip", "", basename(fs))
     cat(paste(file.path(txt, txt, "w001001.adf"), collapse = " "), "\n")

Mosaic them all with GDAL, and do a quick decimation, choose your own percentage. (This is longlat data so we really need more sophisticated resampling to get data evenly distributed on the globe, doable in R but for another time.)

    gdalbuildvrt geoid.vrt n45w180/n45w180/w001001.adf n45w135/n45w135/w001001.adf      n45w90/n45w90/w001001.adf n45w45/n45w45/w001001.adf n45e00/n45e00/w001001.adf n45e45/n45e45/w001001.adf n45e90/n45e90/w001001.adf n45e135/n45e135/w001001.adf n00w180/n00w180/w001001.adf n00w135/n00w135/w001001.adf n00w90/n00w90/w001001.adf n00w45/n00w45/w001001.adf n00e00/n00e00/w001001.adf n00e45/n00e45/w001001.adf n00e90/n00e90/w001001.adf n00e135/n00e135/w001001.adf s45w180/s45w180/w001001.adf s45w135/s45w135/w001001.adf s45w90/s45w90/w001001.adf s45w45/s45w45/w001001.adf s45e00/s45e00/w001001.adf s45e45/s45e45/w001001.adf s45e90/s45e90/w001001.adf s45e135/s45e135/w001001.adf s90w180/s90w180/w001001.adf s90w135/s90w135/w001001.adf s90w90/s90w90/w001001.adf s90w45/s90w45/w001001.adf s90e00/s90e00/w001001.adf s90e45/s90e45/w001001.adf s90e90/s90e90/w001001.adf s90e135/s90e135/w001001.adf 


    gdal_translate geoid.vrt geoid04.tif -outsize 4% 4% -co COMPRESS=LZW -co TILED=YES

Now back to R, read that simplified raster and convert to 3D.

library(rglwidget)
library(raster);library(rgdal)
## Loading required package: sp
## rgdal: version: 1.1-4, (SVN revision 594)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: E:/inst/R/R/library/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: E:/inst/R/R/library/rgdal/proj
##  Linking to sp version: 1.2-1
r <- raster(file.path(gpath, "geoid04.tif"))
extent(r) <- extent(-180, 180, -90, 90)
library(gris)  ## devtools::install_github("mdsumner/gris")  

# Build geometry
geoid0 <- bgl(r, z = r * 5000)
geoid <- geoid0
## convert long-lat-z to XYZ
z <- t(llh2xyz(t(geoid$vb[1:3, ])))


## build map
library(rworldmap)
## ### Welcome to rworldmap ###
## For a short introduction type :   vignette('rworldmap')
data(coastsCoarse)
globe <- gris(coastsCoarse)

## get natural earth image
u <- "http://naturalearth.springercarto.com/ne3_data/8192/textures/2_no_clouds_8k.jpg"
if (!file.exists(basename(u))) download.file(u, basename(u), mode = "wb")
tmp <- setExtent(brick(basename(u)), extent(-180, 180, -90, 90))
projection(tmp) <- "+proj=longlat +ellps=WGS84"
tmp <- projectRaster(tmp, raster(extent(tmp), nrow = nrow(tmp)/8, ncol = ncol(tmp)/8, crs = projection(tmp)))


imname <- gsub("jpg$", "png", basename(u))
if (!file.exists(imname)) writeGDAL(as(tmp, "SpatialGridDataFrame"), imname, drivername = "PNG")

And plot

library(colorspace)
## 
## Attaching package: 'colorspace'
## The following object is masked from 'package:raster':
## 
##     RGB
pal <- colorRampPalette(c("blue", "lightblue", "white", "orange", "red"))(24)
scl <- function(x) (x - min(x, na.rm = TRUE)) / diff(range(x, na.rm = TRUE))
library(rgl)
## 
## Attaching package: 'rgl'
## The following object is masked from 'package:gris':
## 
##     triangulate
shade3d(geoid, col = pal[scl(geoid0$vb[3,geoid$ib]) * (length(pal) - 1) + 1])
# for (i in seq(nrow(globe$o))) {
#     xx <- mkpslg(globe[i, ])
#     xx$XYZ <- llh2xyz(cbind(xx$P, extract(r, xx$P) + 1000 ))
#     lines3d(xx$XYZ[xx$S, ])
# }

subid <- currentSubscene3d()
rglwidget(elementId="geoid")

# imagesphere <- bgl(raster(extent(-180, 180, -90, 90), nrow = 180, ncol = 360))
# tcoords <- xyFromCell(setExtent(tmp, extent(0, 1, 0, 1)), cellFromXY(tmp, t(imagesphere$vb[1:2, ])))
# imagesphere$vb[1:3, ] <- t(llh2xyz(t(imagesphere$vb[1:3, ])))
# 
# shade3d(imagesphere, col = "white", texcoords = tcoords[imagesphere$ib, ], texture = imname)
# 
# subid <- currentSubscene3d()
# rglwidget(elementId="sphere")